SM  Report  86-21 


A  FINITE  ELEMENT  ANALYSIS 
OF  SMALL-SCALE  YIELDING  NEAR 
A  STATIONARY  CRACK  UNDER  PLANE  STRESS 

by 

R.  Narasimhan*  and  A.  J.  Rosakis** 


June  1986 


Division  of  Engineering  and  Applied  Science 
California  Institute  of  Technology 
Pasadena,  CA  91125 


*  Graduate  Research  Assistant. 

**  Assistant  Professor  of  Aeronautics  and  Applied  Mechanics. 


-1- 


ABSTRACT 

A  detailed  finite  element  analysis  of  the  monotonic  loading  of  a  stationary 
crack  is  performed  under  Mode  I  plane  stress,  small-scale  yielding  conditions.  A 
small  strain,  J2  incremental  plasticity  theory  is  employed  and  both  elastic-perfectly 
plastic  and  power  law  hardening  materials  are  considered.  Some  issues  such  as 
the  range  of  dominance  of  the  asymptotic  stress  and  deformation  fields  and  the 
amount  of  non-proportional  loading  near  the  crack  tip,  which  have  received  wide 
attention  in  the  analogous  plane  strain  problem,  are  examined.  Special  attention  is 
devoted  to  the  perfectly  plastic  idealization  by  performing  a  separate  singular  finite 
element  analysis  to  clarify  some  details  about  the  asymptotic  stress  and  deformation 
fields.  The  full-field  numerical  solution  is  used  to  simulate  synthetic  (optical)  caustic 
patterns  at  different  distances  from  the  tip,  which  are  compared  with  experimental 
observations  and  with  asymptotic  analytical  results. 
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1.  INTRODUCTION 

The  stress  intensity  factor  is  a  measure  of  the  intensity  of  the  stress  and  strain 
fields  near  a  crack  tip  in  linear  elastic  fracture  mechanics.  However,  fracture  in 
most  structural  materials,  particularly  low  and  intermediate  strength  metals  is  often 
accompanied  by  plastic  flow  near  the  crack  tip,  invalidating  the  assumptions  of  linear 
elasticity  theory.  Under  certain  circumstances,  the  stress  intensity  factor  can  still 
be  used  to  characterize  the  onset  of  crack  growth,  provided  that  the  plastic  zone 
is  contained  well  within  the  region  of  dominance  of  the  singular  elastic  field.  This 
situation  is  often  referred  to  as  “small-scale  yielding.”  But  when  plastic  flow  takes 
place  over  large  size  scales,  one  is  compelled  to  seek  continuum  solutions  for  crack 
problems  within  the  context  of  an  elastic-plastic  theory. 

HUTCHINSON  (1968a,  1968b)  and  RICE  and  ROSENGREN  (1968) 
performed  the  asymptotic  analysis  for  stress  and  deformation  fields  near  a  mono- 
tonically  loaded  stationary  crack  tip  in  a  power  law  hardening  material  obeying 
a  deformation  plasticity  theory.  The  fact  that  the  value  of  the  J  integral  (RICE, 
(1968a))  provides  a  measure  of  the  intensity  of  the  near-tip  field  in  this  asymptotic 
solution  has  prompted  some  investigators  (e.g.,  BEGLEY  and  LANDES  (1972)) 
to  propose  a  criterion  for  the  onset  of  crack  growth  based  on  the  attainment  of 
a  critical  value  for  J.  This  proposal  has  been  complemented  by  a  wide  range  of 
experimental  data  (e.g.,  LANDES  and  BEGLEY  (1972)). 

In  order  to  characterize  fracture  initiation  based  on  this  single  macroscopic  pa¬ 
rameter,  it  is  imperative  that  the  plastic  singular  fields  of  HUTCHINSON  (1968a, 
1968b)  and  RICE  and  ROSENGREN  (1968)  should  dominate  over  a  length  scale 
that  is  large  as  compared  to  the  fracture  process  zone.  In  this  region,  microstruc- 
tural  processes  such  as  void  nucleation  and  growth,  microcracking,  etc.  take  place. 
The  fracture  process  zone  is  often  believed  to  coincide  with  the  region  near  the  tip, 
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wherein  finite  strain  effects  are  significant.  In  addition  to  the  above  issues,  another 
important  factor  that  has  to  be  considered  is  the  possibility  of  non-proportional 
loading  near  the  tip,  which  would  render  the  deformation  plasticity  theory  (on 
which  the  analysis  of  HUTCHINSON  (1968a,  1968b)  and  RICE  and  ROSENGREN 
(1968)  is  based)  to  be  physically  inappropriate. 

The  above  issues  have  been  examined  by  several  investigators  through  nu¬ 
merical  methods  predominantly  under  the  tensile  plane  strain  mode  of  fracture. 
Accurate  finite  element  studies  with  crack  tip  elements  making  use  of  special  in¬ 
terpolation  functions  to  account  for  the  plastic  strain  singularity  were  conducted 
by  LEVY,  MARCAL,  OSTERGREN  and  RICE  (1971)  and  RICE  and  TRACEY 
(1973)  for  the  perfectly  plastic  case  and  by  TRACEY  (1976)  for  hardening  mate¬ 
rials.  These  studies  modelled  Mode  I  plane  strain,  small-scale  yielding  conditions 
and  employed  an  incremental  plasticity  theory.  They  confirmed  the  validity  of  the 
dominant  fields  of  HUTCHINSON  (1968a,  1968b)  and  RICE  and  ROSENGREN 
(1968)  in  a  region  quite  close  to  the  crack  tip.  McMEEKING  (1977)  performed  a 
finite  element  calculation  to  model  crack  tip  blunting  based  on  a  finite  strain  in¬ 
cremental  plasticity  theory  under  plane  strain,  small-scale  yielding  conditions.  He 
observed  that  finite  strain  effects  become  important  only  for  distances  from  the  tip 
of  the  order  of  2  or  3  times  the  crack  opening  displacement  <5f  (which  will  be  defined 
in  Sec. (4)).  Strong  path  dependence  of  the  J  integral  was  also  noticed  within  this 
region. 

SHIH  and  GERMAN  (1981)  investigated  the  range  of  dominance  of  the  plastic 
singular  fields  for  a  wide  variety  of  specimen  configurations  and  material  proper¬ 
ties  from  contained  yielding  to  fully  plastic  conditions.  They  employed  a  small 
strain  incremental  plasticity  theory  and  confined  their  attention  to  Mode  I  plane 
strain.  McMEEKING  and  PARKS  (1979)  also  investigated  configuration  depen- 


-4- 


dence  within  the  context  of  a  finite  strain  theory  similar  to  (McMEEKING,  (1977)) 
under  large  scale  yielding.  Thus,  substantial  work  under  Mode  I  plane  strain  con¬ 
ditions  has  been  performed  to  provide  a  better  understanding  of  the  mechanics  of 
crack  tip  state  and  also  to  specify  size  requirements  for  specimens  used  in  fracture 
toughness  testing  to  ensure  J  dominance. 

However,  very  little  information  is  available  in  the  literature  pertaining  to 
the  above  issues  under  Mode  I  plane  stress,  despite  its  practical  importance  to 
structural  problems.  A  preliminary  numerical  investigation  was  carried  out  by 
HILTON  and  HUTCHINSON  (1971)  under  plane  stress,  small-(and  large-)  scale 
yielding  conditions  in  which  the  plastic  singular  fields  were  imposed  in  a  small  circle 
near  the  crack  tip.  The  value  of  J  or  some  other  equivalent  plastic  intensity  factor 
was  determined  along  with  the  nodal  displacements  from  the  finite  element  solution. 
SHIH  (1973)  applied  their  method  to  study  combined  Mode  I  and  Mode  II  fracture 
problems  under  both  plane  strain  and  plane  stress.  Both  these  studies  employed 
a  deformation  plasticity  theory  and  considered  power-hardening  materials.  Also, 
the  validity  of  the  asymptotic  solution  of  HUTCHINSON  (1968a,  1968b)  and  RICE 
and  ROSENGREN  (1968)  was  assumed  over  a  length  scale,  which  was  not  known 
a  prion ,  although  this  was  contained  well  within  the  plastic  zone  in  these  numerical 
simulations. 

Some  of  the  issues  mentioned  above,  pertaining  to  the  range  of  dominance  of  the 
asymptotic  fields  and  the  amount  of  non-proportional  loading  near  the  tip,  which 
have  received  considerable  attention  in  the  plane  strain  problem,  have  not  been 
examined  in  plane  stress.  Thus,  detailed  numerical  work  along  the  lines  of  RICE 
and  TRACEY  (1973),  McMEEKING  (1977)  and  SHIH  and  GERMAN  (1981)  is 
required  to  firmly  establish  a  conceptual  understanding  of  fracture  under  plane 
stress  conditions.  This  is  usually  more  complex  than  in  plane-  strain,  primarily 
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because  the  equations  of  plane  stress  plasticity  are  somewhat  more  involved  (e.g., 
HILL  (1983)). 

In  addition  to  the  above  considerations,  a  detailed  numerical  study  of  plane 
stress  fracture  is  important  because  of  the  possibility  of  a  direct  comparison  with 
optical  experimental  methods  such  as  caustics.  This  method,  which  has  been  em¬ 
ployed  to  determine  the  stress  intensity  factor  in  linear  elastic  fracture  problems 
(e.g.,  THEOCARIS  and  GDOUTOS  (1972)),  has  recently  been  extended  to  appli¬ 
cations  in  ductile  fracture  (ROSAKIS,  MA  and  FREUND  (1983);  ROSAKIS  and 
FREUND  (1982)).  A  knowledge  of  the  range  of  dominance  of  the  plastic  singular 
fields  is  of  primary  importance  to  facilitate  a  proper  interpretation  of  experimental 
data  (ZEHNDER,  ROSAKIS  and  NARASIMHAN  (1986)).  Also,  information  from 
full-field  numerical  solutions  would  be  crucial  in  analysing  the  caustics  obtained  in 
regions  outside  the  range  of  dominance  of  any  particular  asymptotic  field. 

In  this  work,  an  elaborate  finite  element  investigation,  with  a  very  fine  mesh 
elucidating  the  details  near  the  crack  tip,  is  undertaken  to  simulate  Mode  I  plane 
stress,  small-scale  yielding  conditions.  No  attempt  has  been  made  in  this  part  of 
the  work  to  incorporate  the  expected  singularity  in  the  strains  by  using  special 
crack  tip  elements.  Computations  have  been  performed  for  materials  obeying  an 
incremental  plasticity  theory  with  no  hardening  and  with  a  power-law  hardening. 
In  Sec. (2),  the  numerical  formulation,  finite  element  scheme,  etc.  are  outlined. 
In  Sec. (3),  stationary  crack  tip  fields  under  plane  stress  are  reviewed.  In  Sec. (4), 
detailed  results  are  presented  for  the  plastic  zones,  stress  and  strain  distributions, 
and  crack  opening  displacement.  Also,  the  path  independence  of  the  J  integral  is 
examined. 

In  Sec. (5),  caustic  patterns  are  simulated  from  the  numerical  solution  at  a  wide 
range  of  distances  from  the  crack  tip  and  are  compared  with  experimental  observa- 
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tions  (ZEHNDER  et  al.  (1986))  and  asymptotic  results  (ROSAKIS  et  al.  (1983)). 
In  Sec. (6),  an  additional  numerical  analysis,  employing  singular  elements  near  the 
crack  tip,  is  performed  for  the  perfectly  plastic  case  in  order  to  examine  the  asymp¬ 
totic  stress  and  deformation  fields.  The  issue  of  sensitivity  of  the  numerical  results 
to  the  near-tip  mesh  design  is  thus  investigated.  It  is  found  that  the  dominant 
strain  field  near  the  tip  for  perfect  plasticity  is  completely  different  from  the  limit 
of  the  singular  solution  of  HUTCHINSON  (1968a,  1968b)  and  RICE  and  ROSEN- 
GREN  (1968)  for  materials  with  low  hardening.  On  the  other  hand,  the  numerical 
results  for  the  near-tip  stress  field  are  in  good  agreement  with  the  slip  line  solution 
of  HUTCHINSON  (1968b).  In  the  light  of  this  observation,  it  is  suggested  that 
the  configuration  dependence  of  crack  tip  deformation  should  be  investigated  under 
plane  stress  in  the  spirit  of  SHIH  and  GERMAN  (1981)  and  McMEEKING  and 
PARKS  (1979).  Such  an  analysis  could  be  complemented  by  experimental  results 
based  on  caustics. 

2.  NUMERICAL  ANALYSIS 

Formulation 

The  Mode  I  plane  stress,  small-scale  yielding  problem  (RICE  (1968b))  was 
modelled  by  considering  a  crack  in  a  domain  R,  which  was  entirely  represented  by 
finite  elements  as  shown  in  Figs,  la  and  b.  Only  the  upper  half-plane  was  considered 
because  of  Mode  I  symmetry.  All  field  quantities  are  referred  to  with  respect  to  an 
orthonormal  frame  {e^eg,^}  centered  at  the  crack  tip.  The  leading  term  in  the 
displacements  of  the  linear  elastic  asymptotic  solution, 

=  ,  (2.1) 

was  specified  as  boundary  conditions  on  the  outermost  boundary  S  of  the 
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Figure  1.  Finite  element  mesh:  a.)  Outer  mesh  b)  Fine  mesh  near  the  crack  tip. 
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domain  1  The  loading  was  applied  through  the  Mode  I  stress  intensity  factor  K[, 
which  occurs  as  an  amplitude  factor  in  Eqn.(2.1). 

The  maximum  extent  of  the  plastic  zone  surrounding  the  crack  tip  was  at  all 
times  within  ^  of  the  radius  of  the  outermost  contour  S,  so  that  the  small-scale 
yielding  condition  was  preserved.  All  plastic  deformation  was  confined  within  the 
active  region  shown  in  Fig.  la,  which  has  a  total  of  1704  four-noded  elements  and 
3549  degrees  of  freedom.  The  large  region  surrounding  this  active  mesh  has  a  total 
of  40  rings  with  56  elements  in  each  ring  and  remained  elastic  throughout  the  entire 
computation.  The  constant  stiffness  of  this  region  was  statically  condensed  using 
a  ring-by-ring  static  condensation  procedure  that  involved  a  partial  forward  Gauss 
reduction  at  each  stage. 

The  cutout  in  Fig.  la  is  a  fine  mesh  region  near  the  crack  tip,  which  is  shown 
in  detail  in  Fig.  16.  This  mesh  was  designed  to  have  small  rectangular  elements 
parallel  to  the  crack  plane  instead  of  being  focussed  at  the  crack  tip.  No  attempt 
has  been  made  to  incorporate  the  singularity  of  the  plastic  strains  by  using  special 
crack  tip  elements  in  this  analysis  (see  Sec. (6)  and  RICE  and  TRACEY  (1973)). 
This  was  because  the  stress  and  strain  fields  at  the  end  of  the  stationary  load 
history  were  used  as  initial  conditions  for  simulating  stable  crack  extension,  which 
will  be  reported  elsewhere.  The  radius  RA  of  the  active  mesh  and  the  radius  of  the 
outermost  boundary  S  are  about  385  times  and  3400  times  the  size  L  of  the  smallest 
element  near  the  crack  tip,  respectively. 

The  Mode  I  symmetry  conditions  that  are  given  by 

<7i2(zi,£2  =  0)  =0  ) 

>,  zi>0,  (2.2) 

u2(xi,a:2  =  0)  =  0  J 

1  Throughout  this  paper,  Greek  subscripts  will  have  the  range  1,2,  while  Latin 
subscripts  will  take  values  1,2,3. 
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were  imposed  by  attaching  stiff  springs  in  the  £2  direction  to  the  nodes  ahead  of 
the  crack  tip.  Traction-free  conditions  were  imposed  on  the  crackflank. 

The  type  of  element  used  was  the  four  noded  isoparametric  quadrilateral,  which 
was  formed  from  four  constant  strain  triangles  with  static  condensation  of  the  in¬ 
ternal  node.  This  element  was  suggested  by  NAGTEGAAL,  PARKS  and  RICE 
(1974)  to  relieve  artificial  mesh-locking  effects  that  occur  under  nearly  incompress¬ 
ible  conditions  in  plane  strain.  However,  this  problem  does  not  arise  in  plane  stress 
because  there  is  a  non-zero  out-of-plane  strain  component  633,  which  is  determined 
in  terms  of  the  in-plane  strain  components  eap. 

Material  Idealization 

The  materials  that  were  numerically  modelled  were  the  elastic-perfectly  plastic 
solids  and  isotropic  power-hardening  solids.  A  small  strain  incremental  plasticity 
theory  was  employed  along  with  the  Huber- Von  Mises  yield  condition  and  the  asso¬ 
ciated  flow  rule.  The  Huber- Von  Mises  yield  condition  for  isotropic  hardening  takes 
the  form, 

ffe.fp)  =  Ffe)  -  a2(ep)  ,  (2.3) 

1/2 

where  F(<r)  =  §S  .  S  and  ep  =  /  (|ep  ep  )  dt  is  the  accumulated  equivalent  plastic 
strain.  In  the  above,  S  is  the  deviatoric  stress  tensor  and  <r(<Tp)  is  defined  by  the 
following  power  hardening  rule: 


eP 


eo 


(2.4) 


For  the  elastic-perfectly  plastic  case,  a  takes  the  constant  value  of  ao,  the  yield 
stress  in  uniaxial  tension.  In  Eqn.(2.4),  eo  is  the  yield  strain  in  uniaxial  tension. 

Within  the  context  of  the  small  strain  flow  theory  of  plasticity,  the  total  strain 
rate  tensor  can  be  decomposed  into  elastic  and  plastic  parts: 


ee  +  ep 


(2.5) 
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The  stress  rate  tensor  &  is  related  to  the  elastic  strain  rate  tensor  ce  through  a 
constant,  isotropic,  positive  definite  elasticity  tensor  C  as, 


a 


(2.6) 


The  plastic  strain  rate  tensor  ep  is  normal  to  the  yield  surface  and  the  flow  rule 
takes  the  form, 

iP  =  jF^=AS,  (2.7) 

where  A  >  0. 

By  using  Eqns.(2.3)-(2.7)  the  constitutive  law  for  material  currently  experienc¬ 
ing  elastic-plastic  deformation  can  be  obtained  as, 


—  ^ijkl^kl  — 


Cijki 


C ijpq  Spq S mnC mnkl 

SrtCrtuvSuv  +  |o2  H 


*kl 


(2.8) 


In  the  above,  H  =  Jp-  and  can  be  obtained  from  (2.4)  for  hardening  solids  and  is 
set  equal  to  zero  for  perfect  plasticity. 

In  the  present  analysis,  (2.3)  and  (2.8)  were  used  along  with  the  plane  stress 
constraint,  which  requires 


°3t  =  0  . 


(2.9) 


By  using  (2.9)  in  (2.8),  an  expression  for  633  can  be  obtained  in  terms  of  iap- 
Finite  Element  Scheme 


A  displacement  based  finite  element  method  was  employed  in  the  analysis.  The 
finite  element  equations  were  derived  from  the  principle  of  virtual  work.  At  a  time 
(t  +  At)  this  takes  the  form, 

I  a(t  +  At) .  <5edA  =  /  T(t  +  At) .  <5uds  .  (2.10) 

Jr  Jd  r 

Here  o(t  +  At)  represents  the  Cauchy  stress  tensor,  which  satisfies  equilibrium  at 
time  (t  +  At)  and  T(t  +  At)  the  imposed  traction  vector  on  the  boundary  dR. 
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Also,  6 u  represents  the  virtual  displacement  vector  that  vanishes  on  the  part  of 
the  boundary  where  the  displacements  are  specified  and  8e  is  the  associated  small 
strain  tensor. 

After  linearizing  about  the  equilibrium  configuration  at  time  t  and  introducing 
the  finite  element  approximation,  the  following  incremental  equilibrium  equations 
are  obtained  in  matrix  form  (e.g.,  BATHE  (1982)): 

KTAU  =  F(t  +  At)  -  P(t)  .  (2.11) 

Here  AU  =  U(t  +  At)  —  U(t)  is  the  vector  of  nodal  point  displacement  increments. 
Also,  Kt  =  fR  BTD  BdA  is  the  tangent  stiffness  matrix  corresponding  to  the  config¬ 
uration  at  time  t,  B,  the  strain  displacement  matrix  (e  =  BU)  and  D,  the  material 
constitutive  matrix.  D  will  be  equal  to  C  for  purely  elastic  response  and  C*  for 
elastic-plastic  material  response.  F(t  +  At)  is  the  vector  of  externally  applied  nodal 
point  loads  at  time  (t  +  At)  and  P(t)  =  fR  BTa(t)dA  is  the  vector  of  nodal  point 
forces  equivalent  to  the  element  stresses  at  time  t. 

In  the  present  analysis,  time  is  only  a  convenient  variable  that  represents  differ¬ 
ent  levels  of  load  intensities.  An  iterative  Newton-Raphson  procedure  (e.g.,  BATHE 
and  CIMENTO  (1980),  BATHE  (1982))  was  employed  in  the  solution  of  the  incre¬ 
mental  equilibrium  Equations  (2.11).  This  method  is  summarized  in  the  Appendix. 

Stress  Computation 

As  was  observed  above,  the  finite  element  scheme  solves  the  displacement  equa¬ 
tions  of  equilibrium  in  an  incremental  fashion.  Hence,  the  constitutive  laws  pre¬ 
sented  earlier  that  deal  with  stress  and  strain  rates  were  used  approximately  to 
relate  small  finite  increments  in  stresses  and  strains.  An  explicit  integration  proce¬ 
dure  also  known  as  the  Tangential  Predictor-Radial  Return  method  was  employed  to 
integrate  the  incremental  stress-strain  law.  As  shown  by  SCHREYER,  KULAK  and 
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KRAMER  (1979),  this  method,  if  used  with  subincrementation  (as  in  the  present 
analysis),  is  very  accurate  for  plane  stress  conditions. 

It  is  important  to  recall  that  the  requirement  of  plane  stress  imposes  a  con¬ 
straint  for  the  out-of-plane  strain  increment  Ae33  in  terms  of  the  in-plane  strain 
increments  Aea/j.  Due  to  this  constraint,  it  is  more  convenient  to  perform  com¬ 
putations  with  stress  and  strain  tensors  instead  of  with  their  deviatoric  parts  as  is 
normally  done  in  plane  strain.  The  method  of  stress  computation  is  outlined  in  the 
Appendix. 

Solution  Strategy 

As  noted  earlier,  the  loading  was  applied  through  the  Mode  I  stress  intensity 
factor  Ki,  which  enters  the  far-held  displacement  boundary  condition  (2.1).  An 
initial  load  step  was  performed  in  which  Ki  was  small  enough  to  ensure  that  all 
the  elements  remained  elastic.  Ki  was  then  scaled  to  cause  incipient  yielding  in  the 
element  nearest  to  the  crack  tip. 

Subsequent  load  steps  were  performed  by  increasing  Ki  by  5-10%  of  the  in¬ 
cipient  value  at  a  time  and  iterating  for  convergence  to  equilibrium.  Each  load 
step  required  typically  3-4  iterations  before  converging  to  an  accepted  equilibrium 
configuration.  Yielding  was  continued  till  the  plastic  zone  surrounding  the  crack 
tip  had  a  maximum  extent  of  about  50  or  100  times  the  smallest  element  size  L  in 
order  to  guarantee  sufficient  resolution  near  the  crack  tip. 

3.  STATIONARY  CRACK  TIP  FIELDS 

Power-Hardening  Solids 

HUTCHINSON  (1968a,  1968b)  and  RICE  and  ROSENGREN  (1968)  investi¬ 
gated  the  asymptotic  stress  and  strain  fields  near  a  monotonically  loaded  stationary 
crack  tip  in  an  elastic-plastic  solid.  The  dominant  singular  term  of  their  analysis 
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will  be  referred  to  as  HRR  in  the  sequel.  In  their  work,  a  J2  deformation  plasticity 
theory  and  a  power-law  hardening  idealization  similar  to  (2.4)  were  assumed. 

The  HRR  analysis  employs  a  small  strain  formulation  and  assumes  a  separable 
form  in  polar  coordinates  r  and  9,  for  the  dominant  term  of  the  solution,  to  obtain, 


v 


J 


1 

n+  1 


«?.  ~ 
*} 


J 


eo 


n+  1 


dij(6,n) 


(3.1) 


_ooeolnr " 

In  (3.1),  (7q  and  eo  are  the  yield  stress  and  strain  in  uniaxial  tension  and  n  is  the 
hardening  exponent.  The  angular  factors  a ty  ( 9 ,  n)  and  e?- ( 9 ,  n )  depend  on  the  mode 
of  loading  and  on  the  hardening  exponent.  The  dimensionless  quantity  In,  which 
is  defined  in  (HUTCHINSON  (1968a)),  decreases  from  5  for  n  =  1  to  about  2.6  for 
n  — ♦  oo  under  plane  stress.  J  in  (3.1)  is  the  value  of  the  J  integral  of  RICE  (1968a). 

For  plane  deformations,  the  J  integral  is  defined  for  any  path  of  integration  T 
by, 

J  =  J  (Wut  -  Vi<TijUjti)ds  ,  (3.2) 


where  W  is  the  local  stress  work  density,  a  unit  vector  normal  to  T  and  vl\  is 
a  particle  displacement  vector.  For  our  purposes,  T  will  denote  an  open  contour 
surrounding  the  crack  tip.  The  integral  (3.2)  has  the  well-known  property  of  path 
independence  for  a  wide  class  of  solids,  including  materials  that  obey  the  deforma¬ 
tion  theory  of  plasticity.  Under  small-scale  yielding  conditions,  J  can  be  evaluvated 
from  contours  taken  in  the  far-field  (K  dominated)  elastic  region  as, 

K2 

J  =  ^  (3.3) 


for  plane  stress.  It  is  important  to  note  that  J  enters  (3.1)  as  an  amplitude  factor 
and  hence  provides  a  unique  measure  for  characterizing  fracture  initiation  at  the 
crack  tip. 
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The  main  limitation  of  the  HRR  analysis  is  the  unknown  range  of  dominance 
(e.g.,  with  respect  to  maximum  extent  of  the  plastic  zone)  of  the  singular  solution. 
This  issue  is  important  since  this  range  of  dominance  should  be  large  as  compared 
with  the  fracture  process  zone  and  the  region  near  the  crack  tip  where  the  small 
strain  plasticity  theory  breaks  down.  From  the  experimental  standpoint,  this  infor¬ 
mation  is  crucial  in  the  proper  interpretation  of  experimental  data  based  on  optical 
measurements  (ZEHNDER  et  al.  (1986)). 

Also,  the  discrepancy  between  the  deformation  theory  and  the  more  appropri¬ 
ate  incremental  theory  of  plasticity  has  to  be  assessed  from  the  context  of  crack  tip 
fields.  In  addition,  another  serious  limitation  that  will  be  pointed  out  later  occurs 
when  the  limit  n  — >  oo  is  taken.  This  is  associated  with  the  change  in  nature  of  the 
governing  equations  in  the  limit  as  the  perfect  plasticity  case  is  approached. 

The  above  issues  will  be  investigated  from  the  point  of  view  of  the  plane  stress 
full-field  numerical  solution  presented  here.  This  solution  simulates  small-scale 
yielding  conditions  and  employs  an  incremental  plasticity  theory. 

Perfectly  Plastic  Solids:  Stress  Field 


For  perfectly  plastic  solids,  the  following  important  assumptions  regarding  the 
asymptotic  nature  of  the  stress  field  are  usually  made, 


dc'l{r'e)  - 

o(l) 


d  e 
do 


d°ij 

dd 
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dr 
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(3.4) 


It  is  important  to  bear  in  mind  that  the  field  equations  for  perfect  plasticity  are 
hyperbolic  2,  while  those  for  hardening  solids  are  elliptic. 

2  For  perfectly  plastic  solids  under  plane  stress,  the  governing  equations  for  the 
stresses  could  be  hyperbolic,  parabolic  or  elliptic  (e.g.,  HILL  (1983)). 
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Equation  (3.4)  can  be  used  to  obtain  asymptotic  forms  of  equilibrium  equations 
and  the  Von  Mises  yield  condition  (RICE  and  TRACEY  (1973)).  These  can  be 
employed  to  show  that  only  two  types  of  asymptotic  plastic  sectors  can  exist  near 
the  crack  tip.  These  are  as  follows  for  plane  stress. 

i)  Centered  Fan  Sector 


In  this  sector,  radial  lines  are  stress  characteristics  and  the  asymptotic  stress 
field  has  the  following  form, 

cr°r(0)  =  rocos(0  -  0O)  ' 

=  2rocos(0  -0O)  >  ,  (3.5) 

°re{9)  =  A>sin(0  -  e0)  , 

where  0o  is  an  arbitrary  constant  angle  and  r0  is  the  yield  stress  in  pure  shear.  • 
it)  Constant  Stress  Sector 

In  this  sector,  the  Cartesian  components  of  the  stresses  are  constant, 


</?(#)  =  ba0  .  (3.6) 

The  constants  bap  are  related  by  the  yield  condition.  Straight  lines  along  which  the 
direct  components  of  the  stress  deviator  S°^  vanish  are  stress  characteristics  (HILL 
(1983)). 

HUTCHINSON  (1968b)  assembled  a  solution  for  the  near-tip  field  comprising 
of  a  combination  of  the  above  sectors  as  shown  in  Fig.  2.  The  region  marked  A 
is  a  centered  fan  sector  extending  from  9  =  0°  to  6  =  79.7°,  while  the  regions  B 
and  C  are  two  constant  stress  sectors,  which  occupy  the  angles  from  9  =  79.7°  to 
9  =  180".  The  stresses  in  Sector  A  are  as  given  by  (3.5)  with  90  =  0.  In  particular, 
it  should  be  noted  that  the  stresses  ahead  of  the  crack  tip  (9  =  0)  are  given  by 


'12 


=  0  . 


°22  —  2fo 


(3.7) 


Figure  2.  Analytical  asymptotic  Reid  near  a  stationary  crack  tip  in  a  perfectly 
plastic  solid  under  plane  stress  represented  by  stress  characteristics.. 
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There  is  also  a  discontinuity  in  the  arr  stress  component  between  the  two  constant 
stress  sectors  B  and  C,  which  is  admissible  as  long  as  the  crack  remains  stationary. 
Perfectly  Plastic  Solids:  Deformation  fields 


As  noted  by  RICE  (1968a)  in  the  case  of  plane  strain,  singularities  in  strains 
result  when  slip  lines  focus  at  a  point  as  in  centered  fan  sectors.  The  displacements 
u;  (or  the  rates  u;  in  a  proper  incremental  formulation  (HILL  (1983)))  are  functions 
of  angle  0  as  the  crack  tip  is  approached  within  centered  fan  sectors  resulting  in 
a  discrete  crack  opening  displacement  at  the  tip.  The  following  assumptions  are 
often  made  (RICE  and  TRACEY  (1973))  about  the  displacements  u;  (or  the  rates 
Ui)  within  centered  fan  sectors, 


Ui(r,0)  ~  u°(0) 

dUi(r’9)  ~um  =  d< 


d& 

dui  /i\ 

r— ~o(i) 


do 


(3.8) 


Since  radial  lines  are  stress  characteristics  in  the  fan,  e£r  is  nonsingular  while 
(or  tee)  and  er<?  (or  Ke)  are  singular  as  0(A)  when  the  crack  tip  is  approached 
within  the  fan.  Thus,  it  is  possible  to  write 

&(*) 


e9d  ~  £o- 


*r6  ~  e0  — 


(0) 


, 


(3.9) 


within  fan  sectors.  The  angular  factors  £^(0)  and  ^g(d)  are  non  unique  and  cannot 
be  determined  from  a  local  analysis.  They  depend  on  a  solution  to  the  entire 
boundary  value  problem.  However,  from  the  flow  rule  , 

lJ  V  2  )  To  ’ 

the  following  relation  can  be  obtained  between  ipde  and  epe, 


(3.10) 


•  p  _  .p  SrQ 

trd  ~  ^69  ~c  > 

•see 
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provided  Sgg  7^  0.  Although  this  equation  strictly  applies  for  the  strain  rates  in  an 
incremental  theory,  it  can  be  used  to  relate  the  total  strains  if  the  stresses  remained 
constant  at  a  material  point  from  the  time  it  was  enveloped  by  the  plastic  zone. 
Hence,  it  is  expected  to  hold  approximately  between  the  asymptotic  angular  strain 
factors  ^pe{9)  and  eprd{6). 

The  dominant  HRR  solution  for  the  stresses  (3.1)  approaches  the  limiting  sli¬ 
pline  distribution  of  perfect  plasticity  as  the  hardening  exponent  n  — *  00.  But  as 
has  been  observed  by  LEVY  et  al.  (1971)  and  RICE  and  TRACEY  (1973)  for 
plane  strain,  one  cannot  in  general  expect  the  HRR  singular  solution  for  the  strains 
as  n  — *•  00  to  be  the  dominant  solution  for  perfect  plasticity  because  of  the  non¬ 
uniqueness  noted  earlier. 

On  the  other  hand,  the  strain  components  are  (in  general)  non-singular  in 
the  constant  stress  sectors  and  the  same  displacement  results  if  the  crack  tip  is 
approached  along  different  radial  lines  in  these  sectors. 

An  expression  for  the  near-tip  J  integral  can  be  obtained  from  the  asymptotic 
form  (3.9)  following  plane  strain  analysis  of  RICE  (1968a).  Taking  the  contour  T 
in  (3.2)  to  be  a  circle  of  radius  r,  one  can  write  (3.2)  as, 


J  =  r  J  | W  cos#  —  a rr  (err  cos  #  —  (er$  —  w)  sin#] 

-  arg  [(er0  +  1 0)  cos  #  —  egg  sin  #]  >d0 


In  the  above  equation,  u>  is  the  rotation,  and 


(3.12) 


u)  =  -€rg  +  o(  -  j,  r  — >  0 


Also, 


e"=o0) 

W  =  Wr  +  o 


r-0, 


(3.13) 


(3.14) 
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where 

r*  (o  \1/2 

Wv  =  J  odep  =  o0ep  «  <70  (-€?•€?•  J 

Taking  r  — *  0  in  (3.12)  and  using  the  asymptotic  equations  (3.5),  (3.9),  (3.13)  and 
(3.14),  one  obtains 

JtiP  =  (^f J0  [2l^re)2  +  i'ePee)2]1/2  cosd  +  iPres'in2d 

+  ep9  sin2  O^dO  , 

where  9*  is  the  maximum  angular  extent  of  the  fan. 

4.  RESULTS  AND  DISCUSSION 

The  computations  were  performed  for  two  levels  of  power  hardening,  n  =  5  and 
9  and  also  for  the  elastic-perfectly  plastic  case,  which  is  referred  to  as  n  =  oo  in  the 
following  discussion  of  the  results.  It  should,  however,  be  noted  that  the  elastic- 
perfectly  plastic  calculation  was  performed  with  H  =  4= j-  =  0  in  the  constitutive 
Equation  (2.8).  The  ratio  of  the  Young’s  modulus  to  the  yield  stress  in  pure  shear 
(E/r0)  was  taken  as  1400  for  the  two  cases  of  power-hardening  and  as  350  for  the 
elastic-perfectly  plastic  calculation.  The  Poisson’s  ratio  was  taken  as  0.3  for  all 
cases. 

Plastic  Zones 

The  plastic  zone  surrounding  the  crack  tip  is  shown  in  Fig.  3  for  the  three 
values  of  hardening  exponent  n.  The  crack  tip  is  situated  at  the  origin  of  the 
coordinate  axes  that  have  been  made  dimensionless  by  the  parameter  (Ki/oo)2. 
This  parameter  has  the  unit  of  length  and  also  contains  a  measure  of  the  far-held 
loading.  Hence,  the  size  of  the  plastic  zone  is  expected  to  scale  with  respect  to  this 
parameter  under  small-scale  yielding  conditions.  A  point  in  the  figure  represents  a 
yielded  integration  station  within  an  element.  It  should  be  noted  that  the  plastic 
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Figure  3.  Plastic  zones  surrounding  the  crack  tip  for  three  levels  of  hardening 
n— 5,9  and  oo. 


Figure  4.  Comparison  of  a )  numerical  with  b)  experimental  plastic  zones  for 
material  with  n  =  9. 
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zone  becomes  less  rounded  and  spreads  more  ahead  of  the  crack  tip  with  decreasing 
hardening  (increasing  n). 

These  plastic  zones  agree  well  in  shape  but  are  slightly  smaller  in  size  as 
compared  with  the  results  of  SHIH  (1973),  who  employed  a  deformation  plas¬ 
ticity  theory  and  used  a  singular  element  near  the  crack  tip.  The  maximum 
extent  of  the  plastic  zone  that  occurs  ahead  of  the  crack  tip  [0  =  0)  is  about 
rp  =  0.22(Ki/<70)2,  0.25(Ki/<7o)2  and  0.29(Ki/<7o)2  for  n  =  5,9  and  oo, 
respectively.  For  comparison,  Shih’s  calculation  indicates  an  rp  of  about 
0.32(Ki/cro)2  for  n  =  25,  and  TADA,  PARIS  and  IRWIN  (1973)  report 
rp  =  ^(Ki/cro)2  for  n  =  oo  based  on  an  approximate  calculation.  The  slightly  larger 
size  of  the  plastic  zone  obtained  by  Shih  could  be  due  to  the  imposition  of  the  HRR 
singular  solution  in  a  small  circle  around  the  crack  tip  in  his  analysis.  The  present 
computation  introduces  no  such  a  prion  constraint. 

In  Fig.  4  the  numerically  obtained  plastic  zone  for  n  =  9  is  compared  with 
the  visual  evidence  of  permanent  plastic  deformation  observed  on  the  surface  of  a 
thin  compact  tension  specimen  (ZEHNDER  et  al.  (1986)).  The  material  used  in 
this  experiment  was  a  4340  carbon  steel  with  a  power-hardening  exponent  of  9  in 
uniaxial  tension.  The  experimental  and  numerical  plastic  zones  agree  well  in  shape 
and  also  in  size  when  the  load  levels  in  the  experiment  were  small  and  there  were 
no  boundary  interaction  effects  (contained  yielding). 

Radial  Distribution  of  Stresses 

The  distribution  of  the  normalized  opening  stress,  022/ To,  along  the  xi  axis 
ahead  of  the  crack  tip  and  within  the  plastic  zone  is  shown  in  Fig.  5.  The  centroidal 
values  of  stress  in  the  row  of  elements  ahead  of  the  the  crack  tip  have  been  used  in 
making  this  plot.  Advantage  has  again  been  taken  of  the  self-similarity  noted  earlier, 
with  the  distance  from  the  crack  tip  being  measured  in  terms  of  the  dimensionless 


Po 


Figure  5.  Radial  distribution  of  opening  stress  ahead  of  the  crack  tip.  The  solid 
lines  represent  the  HRR  asymptotic  distribution. 
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variable  X\/{K\I oq)2 .  The  finite  element  results  agree  to  within  1%  with  the  HRR 
asymptotic  stress  distribution  (3.1),  which  is  shown  by  the  solid  lines  in  the  figure, 
in  the  range  0  <  xx  <  0.08(Ki/ero)2.  For  example,  at  xx  =  0.018(K!/<70)2,  the  ratio 
of  the  finite  element  to  the  HRR  asymptotic  stress  is  3.13/3.14,  2.66/2.67  and 
1.999/2.0  for  n  =  5,9  and  oo,  respectively. 

The  values  given  by  the  HRR  distribution  for  ct22  are  higher  than  their  finite 
element  counterparts  by  about  8%  at  the  elastic-plastic  boundary.  This  is  in  marked 
contrast  to  the  corresponding  result  in  plane  strain  (e.g.,  TRACEY  (1976)),  where 
strong  deviation  of  the  finite  element  solution  from  the  HRR  distribution  was  re¬ 
ported  even  for  small  distances  from  the  crack  tip.  Also,  it  should  be  observed  from 
Fig.  5  that  there  is  only  slight  dependence  of  <r2 2  on  n  for  xx  >  0.15(Ki/cr0)2.  The 
finite  element  values  differ  by  less  than  10%  (with  respect  to  n)  in  this  range. 

The  radial  variation  of  all  the  normalized  stress  components  ahead  of  the  crack 
tip  within  the  plastic  zone  for  the  elastic-perfectly  plastic  case  is  shown  in  Fig.  6. 
The  finite  element  values  near  the  crack  tip  are  in  excellent  agreement  with  the 
asymptotic  slipline  solution  of  Hutchinson  (Fig.  2).  At  xx  =  0.01(Ki/cto)2,  <?n  and 
ct22  are  0.98 r0  and  1.999ro,  respectively,  which  compares  very  closely  with  the  values 
of  t0  and  2r0  given  by  the  slipline  solution  (Eqn.(3.7)).  Also,  Fig.  6  indicates  that 
the  a  11  stress  component  has  a  strong  radial  variation  ahead  of  the  crack  tip,  with 
a  value  at  the  elastic-plastic  boundary  of  about  1.40r0.  This  suggests  curving  of  the 
leading  boundary  of  the  fan  at  moderate  distances  from  the  tip. 

The  plane-stress  Huber- Von  Mises  yield  surface  can  be  represented  by  an  ellipse 
in  principal  stress  space  in  the  following  parametric  form  (HILL  (1983)), 

7T 

o  1  =  2ro  cos(u>  —  — ) 

,  7 r.  1 

0*1  —  2r0  cos  (a;  -f-  —  J  ^  • 

“  =  «  (r,0) 
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Figure  6.  Radial  stress  distribution  ahead  of  the  crack  tip  for  the  perfectly 
plastic  case. 
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Figure  7.  Comparison  of  the  radial  stress  distribution  ahead  of  the  tip  as  given 
by  the  K/  held  ( solid  line)  with  the  finite  element  solution  for  materials  with  a) 
n=5  and  b)  n=9. 
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For  o\  >  cr2,  the  angle  u  varies  in  the  range  0  <  u  <  n.  The  governing  equations 
for  the  stresses  are  hyperbolic  if  f  <  u  <  parabolic  if  u>  =  |  or  and  elliptic 
if  0  <  u;  <  ^  or  ^  <  w  <  7r.  The  value  of  uj(r  — 0,0)  corresponding  to  the 
asymptotic  stresses  (3.7)  is  whereas  the  stresses  at  the  elastic-plastic  boundary 
ahead  of  the  crack  tip  give  w(rp, 0)  «  Thus,  while  the  stress  state  ahead  of 
the  crack  is  parabolic  near  the  tip,  it  appears  to  be  elliptic  at  the  elastic-plastic 
boundary. 

It  is  important  from  the  viewpoint  of  optical  experimental  methods  (such  as 
caustics)  to  determine  the  effect  of  the  crack  tip  plastic  zone  on  the  stress  and 
deformation  fields  in  the  surrounding  elastic  region,  in  order  to  properly  interpret 
the  experimental  data.  To  examine  this  effect,  the  radial  distribution  of  stresses  in 
the  ray  ahead  of  the  crack  tip  is  shown  on  an  expanded  scale  in  Fig.  7  for  the  two 
levels  of  hardening,  n  =  5  and  9.  The  stresses  given  by  the  singular  elastic  solution 
(Ki  field)  are  shown  for  comparison  by  the  solid  line  in  the  figure..  It  is  found  that 
the  C22  stress  component  obtained  from  the  numerical  solution  is  higher  than  that 
given  by  the  singular  elastic  field  at  the  elastic-plastic  boundary  (r  =  rp)  by  more 
than  30%.  However,  the  stress  distribution  undergoes  a  rapid  transition  outside  the 
plastic  zone  and  differs  from  the  Ki  field  by  less  than  8%  for  r  >  1.5rp.  Also,  the 
stress  distribution  in  the  surrounding  elastic  region  seems  to  be  quite  insensitive  to 
the  hardening  level. 

Radial  Distribution  of  Plastic  Strains 

The  radial  variation  of  the  normalized  plastic  strains  e£2/eo  and  £33/60  with 
respect  to  normalized  distance  ahead  of  the  crack  tip  is  shown  in  Fig.  8  for  the  two 
levels  of  power  hardening.  The  HRR  solution  for  the  asymptotic  strain  distribution 
(Eqn.(3.1))  is  shown  by  the  solid  lines  in  the  figure.  The  finite  element  solution, 
although  slightly  smaller  than  the  HRR  distribution  near  the  crack  tip,  appears  to 
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Figure  8.  Radial  variation  of  the  plastic  strains  ahead  of  the  crack  tip  for  ma¬ 
terials  with  a)  n=5  and  b)  n=9  and  comparison  with  HRR  solution  (solid  lines). 
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indicate  the  correct  singular  behaviour  in  the  range  r  <  0.3rp.  It  should  be  recalled 
that  a  very  detailed  mesh  was  used  near  the  crack  tip  (Fig.  16),  and  that  the  plastic 
zone  was  quite  large  as  compared  with  the  smallest  element  size  (at  least  50  times) 
at  the  stage  when  these  results  were  taken.  These  factors  compensate  to  some 
extent  for  the  incorrect  modelling  of  the  singularity  (3.1)  by  our  using  linear  shape 
functions  for  the  crack  tip  elements. 

The  radial  variation  of  the  normalized  plastic  strains  ahead  of  the  crack  tip 
for  the  elastic-perfectly  plastic  case  is  shown  in  Fig.  9.  The  solid  line  in  the  figure 
is  the  limit  of  the  HRR  dominant  singular  solution  for  e^/ eo  for  large  n,  which  is 
given  by  (SHIH  (1973);  ROSAKIS  et  al.  (1983)), 
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The  finite  element  solution  for  the  strains  seems  to  indicate  the  correct  i  variation 
near  the  crack  tip  (r  <  0.04(Ki/<7o)2)  but  is  about  3.3  times  the  values  given  by 


(4.2). 


As  has  already  been  noted  in  Sec. (3),  the  HRR  singular  strain  solution  as 
n  — >  oo,  cannot  (in  general)  be  expected  to  provide  the  dominant  solution  for  per¬ 
fect  plasticity  because  of  the  non-uniqueness  in  strains  associated  with  the  non¬ 
hardening  case.  This  discrepancy  has  also  been  observed  in  Mode  I  plane  strain  by 
LEVY  et  al.  (1971)  and  RICE  and  TRACEY  (1973).  In  this  connection,  it  should 
also  be  mentioned  that  KNOWLES  (1977),  in  working  on  the  finite  anti-plane  shear 
field  near  a  crack  tip  in  an  incompressible  elastic  solid,  with  a  similar  power  law 
behaviour  has  made  an  important  observation.  He  found  that  the  first-  and  second- 
order  terms  in  the  asymptotic  expansion  for  the  displacements  tend  to  become  of 


whi/ 
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Figure  9.  Radial  variation  of  plastic  strains  ahead  of  tip  for  perfectly  plastic 
case.  A  vast  discrepancy  with  the  HRR  singular  solution  for  large  n  (SHIH  (1973); 
ROSAKIS  et  al.  (1983))  is  observed. 
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equal  importance,  as  one  approaches  the  equivalent  of  the  “perfectly  plastic”  case 
in  such  solids.  This  raises  the  question  of  whether  the  limit  as  n  — *  oo  of  the  most 
singular  term  in  the  asymptotic  solution  can  be  considered  separately,  without  ex¬ 
amining  the  limiting  behaviour  of  the  higher-order  terms  of  the  expansion. 

In  order  to  resolve  the  issue  further,  a  separate  finite  element  calculation  for 
the  perfectly  plastic  case  was  performed  under  plane  stress,  small-scale  yielding 
conditions  using  a  focussing  mesh  with  singular  elements  near  the  crack  tip,  similar 
to  the  work  of  RICE  and  TRACEY  (1973).  The  results  of  this  investigation  will  be 
reported  in  Sec. (6).  Finally,  it  should  be  noted  that  the  region  ahead  of  the  crack 
tip,  wherein  the  £  variation  of  the  plastic  strains  was  observed  (r  <  0.04(Ki/cto)2), 
corresponds  to  the  region  of  dominance  of  the  asymptotic  stress  field  (see  Fig.  6). 
Beyond  this  range,  the  front  boundary  of  the  fan  may  tend  to  curve  and  the  £ 
variation  for  the  plastic  strains  may  no  longer  be  valid  (RICE  (1968a,  1968b)). 

Crack  Opening  Displacement 

The  opening  displacement  between  the  crack  faces  as  a  function  of  position 
along  the  crack  flank  is  shown  in  Fig.  10  in  the  nondimensional  form,  <5/(J/oo) 
versus  xi/(Ki/<7o)2,  for  the  three  cases,  n  =  5,9  and  oo.  The  linear  elastic  solution 
corresponding  to  n  =  1  is  also  plotted  for  comparison.  J  in  this  plot  is  the  far-field 
value  given  by  (3.3).  From  the  figure,  it  can  be  observed  that  the  amount  of  blunting 
at  the  crack  tip  increases  with  decreasing  hardening  (or  increasing  n).  There  is  a 
discrete  opening  displacement  at  the  tip  for  the  perfectly  plastic  idealization  because 
of  reasons  stated  in  Sec. (3). 

On  the  other  hand,  the  near-tip  crack  opening  profile  for  the  hardening  cases, 
computed  on  the  basis  of  the  HRR  analysis,  has  the  form  (HUTCHINSON  (1968a); 
RICE  and  ROSENGREN  (1968)), 


6  =  2ii2(r,7r)  ~  (2r(<5t)rl)  "+1 


(4.3) 
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In  this  expression,  St,  which  can  be  written  as 

St  =  —  St(e0,n)  ,  •  (4.4) 

°o 

can  be  approximately  interpreted  as  the  opening  distance  between  the  intercept  of 
two  45°  lines  drawn  back  from  the  crack  tip  to  the  deformed  profile.  This  definition 
was  suggested  by  TRACEY  (1976)  as  a  measure  of  the  crack  tip  displacement  for  a 
hardening  material,  since  6(r  =  0)  =  0  in  this  case,  as  can  be  seen  from  (4.3).  SHIH 
(1981)  has  obtained  the  values  for  St(eo,n)  from  the  HRR  solution  for  both  plane 
stress  and  plane  strain.  It  is  found  (SHIH  (1981))  that  St  is  strongly  dependent  on 
n  and  weakly  on  e<>  Also,  as  n  — »  oo,  <5j  becomes  independent  of  eo  and  takes  the 
value  of  1.0  for  plane  stress. 

From  the  present  finite  element  calculation,  the  value  of  dt/(J/cro)  was  obtained 
by  extrapolating  the  near-tip  crack  profile  to  r  =  0  for  the  non-hardening  case  and  by 
fitting  the  form  (4.3)  to  the  near-tip  profile  for  the  hardening  cases.  SHIH  (1981) 
has  also  computed  the  values  of  <5t/(J/cr0)  for  several  values  of  n  from  his  finite 
element  solution  of  (SHIH  (1973)),  which  as  noted  earlier  employed  a  deformation 
plasticity  theory.  These  results  are  summarized  in  the  following  table. 


Table  1:  Values  of  6tj{J/oo)  for  plane  stress 


a0/E 

n  =  5 

n  =  9 

n  =  25 

n  =  oo 

HRR 

0.0012 

0.40 

0.63 

0.89 

1.0* 

Present  Solution 

0.0012 

0.37 

0.57 

0.85 

SHIH  (1973,  1981) 

0.38 

0.86 

* (extrapolated) 


The  slightly  smaller  values  for  <5t/(J/<7o)  obtained  by  the  present  solution,  as 
compared  to  HRR  for  the  hardening  cases,  can  be  accounted  partially  by  some 
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discrepancy  between  flow  theory  and  deformation  theory  as  explained  below.  But 
the  difference  between  the  present  perfect  plasticity  calculation  and  the  HRR  non¬ 
hardening  limit  is  because  the  latter  is  unable  to  provide  complete  information 
regarding  the  most  singular  term  for  the  strains  in  the  asymptotic  solution  for 
perfect  plasticity,  as  described  above.  This  discrepancy  has  also  been  observed  in 
plane  strain.  The  published  numerical  results  (e.g  ,  SHIH  (1981))  for  <5t / ( J / cr0 ) 
under  plane  strain,  small-scale  yielding  conditions  for  the  perfectly  plastic  case 
range  from  0.63-0.66,  whereas  the  HRR  non-hardening  limit  is  0.78. 

J  integral  calculations 

In  order  to  assess  the  difference  between  the  present  incremental  formulation 
and  the  deformation  plasticity  theory,  the  path  independence  of  the  J  integral  was 
checked.  The  J  integral  (3.2)  was  computed  for  the  hardening  materials  along 
several  contours  surrounding  the  crack  tip,  which  passed  through  the  centroids 
of  the  elements.  The  near-tip  contours  enclosing  the  crack  tip  were  rectangular, 
while  the  far-field  contours  were  circular,  in  keeping  with  the  structure  of  the  mesh 
(Fig.  1).  The  integrand  in  (3  2)  was  calculated,  using  the  averaged  values  of  stresses 
and  strains  at  the  centroids  of  the  elements  lying  in  the  contour  path,  and  the 
integration  was  carried  out  numerically  using  Gauss  quadrature.  It  was  found  that 
very  near  the  crack  tip  (r  <  0.04rp)  there  was  a  small  amount  of  path  dependence. 
However,  after  some  distance  away  from  the  crack  tip,  the  calculated  J  value  was 
virtually  indistinguishable  from  the  remotely  applied  value  (3.3). 

For  a  contour  with  an  average  radius  r  =  0.012(Ki/cro)2,  the  ratio  of  the  cal¬ 
culated  J  value  to  the  remotely  applied  J  was  0.96  and  0.95  for  n  =  5  and  9,  re¬ 
spectively.  For  contours  with  average  radius  r  >  0.05(Ki/<7o)2,  the  calculated  J 
value  was  smaller  than  the  applied  J  by  less  than  1%.  While  the  5%  difference 
for  the  near-tip  contours  is  within  the  realm  of  errors  in  the  discretization  proce- 
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dure  and  in  the  numerical  integration  of  (3.2),  it  also  suggests  small  amounts  of 
non-proportional  loading  experienced  by  a  material  particle  from  the  time  it  was 
enveloped  by  the  plastic  zone.  For  the  elastic-perfectly  plastic  material,  our  accu¬ 
rate  numerical  solution  of  Sec. (6)  was  used  to  estimate  the  near-tip  J  integral,  and 
its  discussion  will  be  deferred  till  then. 

In  order  to  further  check  for  discrepancy  between  the  two  plasticity  theories, 
e^/eo  was  calculated  for  the  hardening  materials  at  the  centroids  in  the  row  of  ele¬ 
ments  ahead  of  the  crack  tip  by  substituting  the  averaged  stresses  in  these  elements 
into  the  expression  given  by  the  J2  deformation  theory.  The  plastic  strain  given 
by  the  J2  deformation  theory  was  about  5%  higher  at  r  =  0.012(Ki/cro)2  than  the 
corresponding  value  given  by  the  incremental  formulation  that  was  reported  earlier 
(Fig.  8).  This  difference  progressively  diminished  as  the  distance  from  the  crack  tip 
increased,  and  it  was  less  than  1%  for  r  >  0.1(Ki/<7o)2. 

5.  NUMERICAL  SIMULATION  OF  CAUSTICS 

Introduction 

The  optical  experimental  method  of  caustics  has  been  applied  to  the  study  of 
linear  elastic  fracture  problems  and  to  the  direct  measurement  of  the  stress  intensity 
factors  (e.g.,  THEOCARIS  and  GDOUTOS  (1972);  ROSAKIS  and  ZEHNDER 
(1985)).  This  method  was  recently  extended  to  the  measurement  of  the  J  integral 
in  ductile  fracture  (ROSAKIS  et  al.  (1983);  ROSAKIS  and  FREUND  (1982))  on 
the  basis  of  the  validity  of  the  plane  stress,  HRR  asymptotic  solution. 

Under  conditions  of  small-scale  yielding,  the  singular  elastic  field  dominates 
well  outside  the  plastic  zone.  Inside  the  plastic  zone,  very  near  the  crack  tip,  the 
HRR  field  dominates.  In  the  transition  region  between  these  two  fields,  no  analyti¬ 
cal  solution  is  available.  This  limits  the  applicability  of  caustics,  and  the  conditions 
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under  which  the  results  reported  in  (ROSAKIS  et  al.  (1983)  and  ROSAKIS  and 
FREUND  (1982))  are  valid  are  uncertain.  Also,  errors  may  be  caused  in  the  mea¬ 
surement  of  Kj  based  on  the  caustics  obtained  from  the  elastic  region  surrounding 
the  plastic  zone.  This  is  because  the  crack  tip  plastic  zone  affects  the  caustic  pat¬ 
terns,  and  an  analysis  based  on  the  Kj  field  may  be  erroneous. 

In  this  section,  the  full-field  numerical  solution  under  small-scale  yielding  is 
used  to  generate  simulated  caustic  patterns.  The  numerical  caustics  are  com¬ 
pared  with  the  corresponding  patterns  observed  from  experiments  (ZEHNDER  et 
al.  (1986)).  The  analysis  of  caustics  based  on  the  numerical  results  is  not  limited 
by  the  assumption  of  the  validity  of  any  particular  asymptotic  field.  Finally,  quali¬ 
tative  and  quantitative  comparisons  of  the  simulated  caustics,  obtained  at  various 
distances  from  the  crack  tip,  are  made  with  the  corresponding  results  based  on  the 
near-tip  HRR  analysis  and  the  remotely  applied  K[  field. 

The  Method  of  Caustics 

Consider  a  set  of  parallel  light  rays  normally  incident  on  a  planar,  reflective 
specimen  that  has  been  deformed  by  tensile  loading.  Due  to  the  deformed  shape 
of  the  specimen,  an  envelope  in  space  called  the  “caustic  surface”  is  formed  by 
the  virtual  extension  of  the  reflected  light  rays  (Fig.  11).  The  intersection  of  this 
surface  with  a  plane  located  at  a  distance  zo  behind  the  specimen  is  called  the 
“caustic  curve”  and  it  bounds  a  dark  region  called  the  “shadow  spot.” 

Let  (xj , X2)  be  a  coordinate  system  on  the  specimen  surface  centered  at  the 
crack  tip  and  (Xi,X2),  a  system  translated  by  a  distance  zo  behind  the  specimen 
surface.  Then  the  mapping  of  a  point  (xi,X2)  on  the  specimen  surface  to  a  point 
(Xi,X2)  on  the  plane  at  zo  due  to  reflection  of  a  light  ray  may  be  described  by 
(ROSAKIS  and  ZEHNDER  (1985)), 


xQ  +  2z0 


dU3(xi,X2) 

d\a 


(5.1) 
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Figure  11.  Formation  of  caustic  due  to  reflection  of  light  from  a  polished,  de¬ 
formed  specimen  surface. 
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The  locus  of  points  on  the  specimen  surface  at  which  the  Jacobian  determinant  of 
the  mapping  (5.1)  vanishes  is  called  the  “initial  curve.”  While  points  on  the  initial 
curve  map  onto  the  caustic  curve,  all  points  both  inside  and  outside  the  initial  curve 
map  outside  the  caustic.  The  position  of  the  initial  curve  may  be  varied  by  changing 


z0. 


For  a  stationary  crack  under  small-scale  yielding  conditions,  if  the  initial  curve 
is  chosen  to  fall  well  outside  the  plastic  zone  and  within  the  region  of  validity  of  the 
Ki  field  (large  values  of  zo),  then  the  resulting  caustic  curve  will  be  an  epicycloid 
(Fig.  12a).  In  such  a  case,  Ki  is  related  to  the  caustic  diameter  D  (which  is  the 
maximum  width  of  the  caustic  in  the  X2  direction)  by  (ROSAKIS  and  ZEHNDER 
(1985)), 


Kt  = 


ED5/2 


(5.2) 


10.7z0t'h  ’ 

where  h  is  the  specimen  thickness.  The  initial  curve  is  circular  and  its  radius  r0  is 
given  by 


r0  =  0  316D  . 


(5.3) 


On  the  other  hand,  if  the  initial  curve  is  chosen  to  fall  well  inside  the  plastic 
zone  and  within  the  region  of  dominance  of  the  HRR  field  (very  small  values  of 
Zo),  then  its  shape  as  deduced  by  ROSAKIS  et  al.  (1983)  will  no  longer  be  circular. 
In  such  a  case,  the  radius  ro  of  the  point  on  the  initial  curve  that  maps  to  the 
maximum  value  of  X2  on  the  caustic  curve  is  given  by 


r0  =  0.385D 


(5.4) 


for  a  hardening  exponent  n  of  9.  Also,  the  value  of  the  J  integral  may  be  obtained 
from  the  caustic  diameter  D  as  (ROSAKIS  et  al  (1983)), 
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Figure  12.  Predicted  caustic  shapes  based  ona)  Kr  Held  and  b)  HRR  asymptotic 
Reid  for  n=9. 
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where  Sn  is  a  numerical  factor  dependent  on  n.  Caustic  curves  thus  obtained  from 
the  HRR  field  for  several  values  of  the  hardening  exponent  are  given  by  ROSAKIS 
et  al.  (1983).  A  typical  caustic  for  n=9  is  shown  in  Fig.  126. 

Results  and  Discussion 

The  discrete  values  of  the  out-of-plane  displacement  U3  obtained  from  the  nu¬ 
merical  solution  at  the  centroids  of  the  elements  were  smoothed  using  a  least-squares 
finite  element  scheme  as  advocated  by  HINTON  and  CAMPBELL  (1974).  The  sur¬ 
face  thus  generated  is  shown  in  Fig.  13  for  a  material  with  a  hardening  exponent  of 
9.  Caustic  patterns  were  simulated  by  mapping  light  rays  point  by  point  from  this 
smoothed  surface  using  Equation  (5.1)  for  different  values  of  zo- 

The  sequence  of  caustics  simulated  from  the  finite  element  solution  for  different 
values  of  zo  is  shown  in  Fig.  14  for  a  material  with  n  =  9.  The  parameter  in 
the  figure  is  the  ratio  of  the  initial  curve  size  to  the  maximum  plastic  zone  extent. 
The  initial  curve  size  ro  was  estimated  approximately  by  using  Equation  (5.4)  for 
caustics  from  within  the  plastic  zone  and  by  Equation  (5.3)  for  caustics  from  outside 
the  plastic  zone.  It  is  seen  from  the  figure  that  for  ^  =  0.19,  the  simulated  caustic 
agrees  in  shape  with  the  caustic  predicted  by  the  HRR  field,  which  is  shown  in 
Fig.  126.  When  =  1.3,  the  numerically  simulated  caustic,  Fig.  14/,  agrees  with 

rp 

the  caustic  predicted  using  the  elastic,  Ki  field  (Fig.  12a). 

A  sequence  of  photographs  of  caustics  (ZEHNDER  et  al.  (1986))  obtained  from 
the  tensile  loading  of  a  thin  compact  tension  specimen  of  4340  carbon  steel  is  shown 
in  Fig.  15.  The  experimental  details,  specimen  dimensions,  etc.  are  described  by 
ZEHNDER  et  al.  (1986).  On  comparing  Figs.  14  and  15  it  is  seen  that  in  both 
cases  there  is  a  transition  from  an  “HRR  caustic”  to  an  “elastic  caustic”  as  ^  goes 
from  0.19  to  1.4.  The  transition  away  from  the  HRR  caustic  appears  to  take  place 
slightly  sooner  in  the  numerical  model  (around  ^  =  0.3)  than  in  the  experiment 

ri> 


Figure  13.  Smoothed  out-of-plane  displacement  field  for  n=9. 
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Figure  14.  Sequence  of  caustics  simulated  from  the  numerical  solution. 
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(ZEHNDER  et  al.  (1986)). 
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(around  j-jj-  =  0.35).  However,  the  general  trend  is  similar  in  both  cases. 

It  is  found  that  both  the  numerical  and  experimental  caustics  retain  the  shape 
predicted  by  the  Ki  field  even  for  p2-  as  small  as  1.0.  Thus,  the  effect  of  the  plastic 
zone  cannot  be  judged  by  mere  observation  of  the  caustic  shape.  The  reason  for  the 
invariance  in  shape  of  the  caustics  is  explained  by  examining  the  angular  variation 
of  the  sum  (an  +  022),  °f  the  direct  stress  components  (as  given  by  the  numerical 
solution),  at  different  distances  outside  the  plastic  zone  as  shown  in  Fig.  16.  It  is 
seen  that  the  sum  (an  +  a22)  generally  follows  the  angular  distribution  given  by 
the  Ki  field,  which  is  shown  by  the  solid  line  in  the  figure  even  for  ^  as  small  as 

rP 

1.2.  However,  the  individual  stress  components  show  more  deviation  from  those  of 
the  Ki  field  for  small  values  of  t2-.  This  observation  is  important,  since  the  caustic 
shape  depends  on  the  angular  variation  of  the  out-of-plane  displacement  component 
U3,  which  in  the  elastic  region,  is  proportional  to  (an  +  022)  under  plane  stress. 
Thus,  it  is  not  surprising  that  the  caustic  shape  resembles  the  “elastic  caustic”  for 
^  as  small  as  1.0. 

rP 

The  numerical  caustics  were  simulated  for  a  fixed  value  of  Ki  (or  the  far- 
field  value  of  J  as  given  by  (3.3))  by  varying  z0  in  the  optical  mapping  Equation 
(5.1).  The  relationship  between  the  diameter  D  of  the  simulated  caustics  and  the 
remotely  applied  J  value  is  shown  in  non-dimensional  form  in  Fig.  17.  The  inverse 
of  the  abscissa  in  the  figure  is  an  indication  of  the  initial  curve  size  or  the  distance 
from  the  crack  tip  at  which  the  information  about  the  deformation  field  is  being 
scrutinized.  Thus  a  very  small  abscissa  value,  (large  Zo  or  small  J)  implies  that  the 
initial  curve  is  far  away  from  the  tip.  A  very  large  abscissa  value,  on  the  other  hand, 
implies  that  the  curve  is  very  near  the  tip,  probably  within  the  range  of  dominance 
of  the  HRR  field.  The  bars  on  the  numerical  results  indicate  the  uncertainty  in 
determining  the  initial  curve  due  to  discretization  of  the  finite  elements. 
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Figure  16.  Angular  distribution  of(on  +<722)  for  different  distances  from  the 
tip.  The  solid  line  is  the  distribution  given  by  the  Ki  field. 


{  Numerical 

-  Linear  Elastic 

- Asymptotic,  HRR 


Figure  17.  Relationship  between  caustic  diameter  and  the  J  integral  as  obtained 


from  the  numerically  simulated  caustics. 
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The  solid  line  in  the  figure  represents  the  variation  of  caustic  size  in  the  Ki 
dominated  region  as  given  by  (5.2)  with  u  —  0.3.  The  dashed  line  gives  the  rela¬ 
tionship  for  the  caustics  from  the  HRR-dominated  region  (5.5).  As  can  be  observed 
from  this  figure,  the  numerical  results  approach  the  elastic  relation  (5.2)  for  small 
abscissa  values  and  the  relation  (5.5)  obtained  from  the  HRR  solution  for  large 
values  of  the  abscissa.  In  the  intermediate  region  there  is  a  transition  from  one 
distribution  to  the  other. 

6.  SINGULAR  FINITE  ELEMENT  ANALYSIS 

Introduction 

In  this  section,  a  detailed  investigation  of  the  perfectly  plastic  case  will  be 
presented,  with  the  view  of  examining  closely  the  discrepancy  between  the  numerical 
results  for  the  near-tip  strains  and  the  corresponding  term  of  the  HRR  solution 
(non-hardening  limit),  which  was  noted  in  Sec. (4).  For  this  purpose,  a  singular 
finite  element  analysis  similar  to  the  plane  strain  work  of  RICE  and  TRACEY 
(1973)  was  carried  out  under  Mode  I  plane  stress,  small-scale  yielding  conditions. 
A  ring  of  focussed  isosceles  triangular-shaped  elements  was  used  near  the  crack  tip 
in  this  computation.  This  mesh  design  is  different  from  the  fine  mesh  employed  in 
the  earlier  analysis  (Fig.  16).  Thus,  the  issue  of  sensitivity  of  the  numerical  results 
presented  earlier  in  Sec. (4)  to  the  near-tip  mesh  design  was  also  examined  through 
this  section  of  our  work. 

Numerical  Scheme 

The  near-tip  elements  that  were  employed  here  provide  a  capability  for  non¬ 
uniqueness  of  displacement  at  the  crack  tip  (LEVY  et  al.  (1971);  RICE  and 
TRACEY  (1973)),  which  is  the  fundamental  feature  of  the  i  plastic  strain  sin¬ 
gularity  within  centered  fan  regions  (Sec. (3)).  This  was  achieved  by  treating  the 
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triangular  elements  at  the  crack  tip  as  degenerate  isosceles  trapezoids  that  have 
a  total  of  four  nodes  (one  at  each  vertex)  with  two  nodes  coinciding  at  the  crack 
tip  (Fig.  18).  The  coincident  nodes  at  the  crack  tip  were  constrained  to  move  as 
a  single  point  till  the  load  level  at  which  incipient  yielding  was  detected  in  one  of 
the  near-tip  elements.  A  special  shape  function  (RICE  and  TRACEY  (1973))  was 
used  up  to  this  load  level  to  provide  the  crack  tip  elements  the  capability  to  model 
the  dominant  elastic  strain  singularity.  Subsequently,  the  coincident  nodes  were 
allowed  to  move  independently  and  the  crack  tip  elements  modelled  the  j  plastic 
strain  singularity. 

The  mapping  of  a  four-noded  rectangle  to  a  triangle  (Fig.  18)  can  be  described 


by 


x  _  xi(j  ~  jKj  ±Jl)  +  xj  (1  ~  Q(1  ~  rj)  +  „k(l  +  Q(1  -  y) 
4  ~  4  ~  4 


T  x 


(1  +  C)(i  +  n) 


(6.1) 


with  the  constraint  x1  =  xj.  Here  (£,77)  is  the  natural  coordinate  system  for  the 

element  and  (xl5X2)  is  a  global  coordinate  system  centered  at  the  crack  tip.  The 

inverse  mapping  of  (£,77)  in  terms  of  a  local  Cartesian  coordinate  system  (s,t),  and 

a  local  polar  coordinate  system  (r,t/>)  for  the  element  is  given  by  (Fig.  18), 

2s  4 

-  1 

(6.2) 


6  = 


so 


77  = 


t/s 


tan  rp 


(to/so)  tan  a  ) 

The  elastic  singularity  element  has  the  shape  function  (RICE  and  TRACEY  (1973)), 


+  uK 


u  -  */) 

2 


+  u 


(I  +  77) 


(6.3) 


Here  u1J  represents  the  unique  displacement  of  the  crack  tip  nodes  i  and  j.  The  above 
element  correctly  models  the  -y/r  variation  in  the  leading  term  for  the  displacements 
of  the  linear  elastic  solution.  Also,  displacement  compatibility  is  satisfied  along  the 
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Figure 


18.  Typical  near-tip  element  used  in  the  singular  finite  element  analysis. 
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edges  i-1  and  j-k  (77  =  ±1)  with  the  adjacant  singular  elements  and  along  the  edge 
1-k  (£  =  1)  with  the  conventional  four-noded  isoparametric  element  that  is  joined 
there. 

As  was  first  pointed  out  by  LEVY  et  al.  (1971),  the  mapping  of  any  four- 
noded  isoparametric  element  to  a  triangle  leads  to  a  L  strain  variation  provided 
that  the  coincident  nodes  are  permitted  to  have  different  displacements.  The  crack 
tip  displacement  for  such  an  element  is  given  by  (Fig.  18), 


u(-1,t?) 


(ul+uj) 

2 


(6.4) 


Following  the  notation  of  (3.9)  and  neglecting  the  elastic  strains  that  are  bounded, 
it  can  be  shown  from  (6.4)  that 


fP  - 
e 66  ~ 


er0  — 


E  (  sec  ip 


oo  \  2  tan  a 
E  (  sec  ip 


Cq  \  4  tan  a 


[— (us  —  uJJ  tan  ip  +  (uj.  —  uJt)] 
[«  ~  ut)  +  K  ~  ut)  t&nrp] 


(6.5) 


where  us  and  ut  are  the  displacement  components  in  the  local  (s,t)  Cartesian  co¬ 
ordinate  system  and  ip  is  the  angle  measured  in  the  local  (r,  ip)  polar  coordinate 
system  (Fig.  18)  for  the  element.  It  should  be  noted  that  the  right-hand  side  of 
(6.5)  is  a  first-order  finite  difference  approximation  to  ep d(ip )  and  epe(ip).  Also,  it 
should  be  noted  that  if  the  two  coincident  nodes  displace  as  a  single  point,  so  that 
u1  =  uJ,  then  this  element  behaves  as  an  ordinary  constant  strain  triangle. 

The  mesh  employed  in  this  analysis  was  similar  to  the  one  used  by  LEVY  et 
al.  (1971).  Only  the  upper  half-plane  was  considered  because  of  symmetry.  The  ac¬ 
tive  mesh  consisted  of  20  rings  with  radii  of  L,  (1.5)2L,  (2.0)2L, . . . ,  (9.5)2L,  (10.0)2L 
and  115L.  These  were  divided  by  25  rays  at  equal  angular  intervals  of  7.5°,  giving 
a  total  of  525  nodes  (including  25  coincident  crack  tip  nodes)  and  480  elements  in 
the  active  mesh.  The  region  outside  consisted  of  14  rings  with  24  elements  in  each 


-51- 


ring  and  always  remained  elastic.  Static  condensation  was  employed  in  this  region 
as  described  in  Sec. (2).  The  radius  of  the  outermost  boundary  S  on  which  the  dis¬ 
placement  boundary  condition  (2.1)  was  specified  was  645L.  The  loading  process 
was  stopped  when  the  maximum  plastic  zone  extent  was  about  of  the  radius  of 
the  outermost  boundary  S,  so  that  the  small-scale  yielding  condition  was  preserved. 
The  symmetry  condition  (2.2)  on  the  9  =  0  ray  and  the  traction-free  condition  on 
the  9  —  7r  ray  were  enforced. 

Every  near-tip  element  was  composed  of  three  subelements  (RICE  and 
TRACEY  (1973)),  each  extending  to  one-third  of  the  height  of  the  element.  A  nine- 
point  numerical  integration  scheme  was  employed  to  integrate  the  element  stiffness 
matrix,  with  integration  stations  at  (£,??  =  -|,0,  §)  and  weighting  factors  of  | 
of  the  area  of  the  element.  For  the  isoparametric  elements  outside  the  innermost 
ring,  the  two-by-two  Gauss  quadrature  scheme  was  used.  The  solution  strategy  was 
the  same  as  that  described  in  Sec.  (2)  with  the  additional  modifications  mentioned 
earlier  in  this  section. 

Results  and  Discussion 

It  can  be  shown  by  substituting  the  dominant  term  of  the  elastic  solution  for 
the  stresses  into  the  plane  stress  Von  Mises  yield  condition  that  incipient  yielding 
will  occur  at  an  angle  of  arccos(^)  «  70.5°.  Also,  the  value  of  the  load  parameter, 
Kjj7(cr0 ^/27rry),  calculated  from  the  analytical  solution  is  0.866  for  initial  yielding  at 
a  radius  of  ry.  Incipient  yielding  occurred  in  the  present  finite  element  computation 
in  the  subelement  between  67.5°  and  75°  with  a  mean  angle  of  71.25°.  The  value  of 
Kf/(o,o y/2irry)  was  0.83,  which  is  in  good  agreement  with  the  analytical  prediction. 

The  radial  distribution  of  stresses  along  the  ray  ahead  of  the  crack  tip  at  incip¬ 
ient  yield  is  shown  in  Fig.  19  in  the  nondimensional  form,  o/tq  versus  r/(Ki/<7o)2. 
The  stresses  given  by  the  finite  element  solution  are  in  excellent  agreement  with  the 
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Figure  19.  Radial  distribution  of  stresses  ahead  of  crack  tip  at  incipient  yielding. 
Solid  line  is  the  singular  elastic  solution. 
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dorainant  elastic  solution,  which  is  shown  by  the  solid  line  in  the  figure.  Also,  the 
angular  distribution  of  stresses  within  the  crack  tip  elements  compared  closely  with 
the  analytical  solution. 

The  plastic  zone  at  the  end  of  the  stationary  load  history  is  shown  in  nondi- 
mensional  coordinates  in  Fig.  20.  This  compares  very  well,  in  overall  features,  with 
the  plastic  zone  obtained  in  the  earlier  analysis  (Fig.  3).  The  maximum  plastic  zone 
extent  is  about  rp  =  O.28(Ki/0o)2  ahead  of  the  crack  tip.  In  the  sub-elements  near¬ 
est  to  the  crack  tip,  yielding  spread  only  from  9  =  0  to  75°,  which  is  in  approximate 
agreement  with  the  centered  fan  region  of  Fig.  2. 

The  radial  stress  distribution  ahead  of  the  crack  tip  within  the  plastic  zone 
also  appeared  similar  to  the  variation  reported  earlier  in  Fig.  6.  In  the  subelement 
nearest  to  the  crack  tip  that  occupies  the  angular  range  from  9  =  0  to  7.5°,  the 
stresses  crn  and  022  reached  the  constant  values  0.99ro  and  1.999ro,  respectively, 
which  agrees  very  well  with  the  analytical  asymptotic  limit  (3.7).  Once  again,  a 
strong  radial  variation  in  the  01 1  stress  component  was  observed  along  the  9  =  0 
ray,  with  a  value  at  the  elastic-plastic  boundary  of  1.40ro. 

The  angular  distribution  of  the  normalized  stress  component  crgg/ro,  within  the 
subelements  nearest  to  the  tip,  is  shown  in  Fig.  21  along  with  the  slip  line  solution 
(solid  line)  of  HUTCHINSON  (1968b).  The  finite  element  solution  shows  good 
agreement  with  the  analytical  distribution  in  the  angular  range  0  <  9  <  80°,  which 
corresponds  to  the  centered  fan  region  in  Fig.  2.  This  was  typical  of  the  other  two 
stress  components  oTg  and  orr  as  well,  with  orr  showing  more  deviation  from  the. 
analytical  solution  a s  9  —>■  80°.  This  result  is  consistent  with  the  fact  that  the  two 
constant  stress  sectors  in  Fig.  2  were  not  detected  by  the  finite  element  solution. 
Also,  the  numerical  result  suggests  that  within  the  fan,  the  focussing  of  the  slip 
lines  may  occur  very  close  to  the  crack  tip  in  the  angular  range  65°  <  9  <  80°. 
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Figure  20.  Plastic  zone  for  the  perfectly  plastic  case  obtained  from  the  singular 
finite  element  analysis. 
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Figure  21.  Angular  distribution  of  crgg/ro  within  sub-elements  closest  to  the  crack 
tip  at  the  end  of  the  loading  process.  Solid  line  is  the  analytical,  asymptotic  solution 
of  HUTCHINSON  (1968b). 
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The  normalized  crack  tip  opening  displacement  <St/(J/<To)>  where  J  is  the  re¬ 
motely  applied  value  of  the  J  integral,  was  calculated  based  on  the  crack  tip  node 
lying  on  the  9  =  n  ray.  It  increased  from  zero  at  incipient  yield  (Kj  =  Kj)  to  a  con¬ 
stant  value  of  0.84  at  Ki  «  3.5Kp  This  value  did  not  change  during  the  subsequent 
part  of  the  loading  process.  The  variation  in  <5t/(J/<70)  during  the  initial  phase 
of  the  loading  process  occurred  since  the  plastic  zone  was  not  fully  developed.  It 
should  be  noted  that  this  quantity  is  in  excellent  agreement  with  the  value  reported 
in  Table  1,  which  was  calculated  on  the  basis  of  the  earlier  analysis. 

The  displacements  of  the  crack  tip  nodes  were  substituted  into  Eqn.(6.5),  with- 
-0  =  0  (corresponding  the  mean  angle  of  the  near-tip  element),  to  determine  the 
angular  factors  e00(9)  and  ^0{9)  of  the  dominant  i  strain  singularity  (3.9).  In 
order  to  compare  with  the  dimensionless  angular  factors  ^j{9,n)  given  by  the  HRR 
analysis  (Eqn.(3.1))  for  large  n,  the  functions  e00(9)  and  e(^(0)  obtained  from  the 
present  finite  element  calculation  for  the  perfectly  plastic  case  were  normalized  as 
follows, 


4(«)  =  7^71  ^ 


~v 


W)  = 


{KijcoY 


:In 


(6.6) 


(Kj/aoY 

Here  In  is  taken  as  2.6  corresponding  to  n  — >  oo  in  the  HRR  solution.  The  functions 
thus  obtained  are  shown  along  with  the  HRR  distribution  for  n=25  (which  is  given 
by  SHIH  (1973))  in  Fig.  22.  It  can  be  seen  that  the  two  angular  functions  are 
completely  different.  It  is  interesting  to  note  that  the  numerical  solution  for  the 
perfectly  plastic  case  under  small-scale  yielding  conditions  gives  vanishingly  small 
values  for  the  angular  factors  of  the  dominant  i  strain  singularity  for  9  >  45°, 
although  the  slip  line  solution  of  Fig.  2  shows  a  centered  fan  extending  from  9  =  0 
to  about  80°. 

It  is  found  that  the  angular  factors  lp00  and  e^,  obtained  from  the  numerical 
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Figure  22.  The  angular  factors  of  the  1/r  plastic  strain  singularity  obtained  from 
the  numerical  solution  on  the  basis  of  the  non-unique  crack  tip  displacement.  Solid 
line  is  the  variation  given  by  the  HRR  solution  for  n=25  (SHIH  (1973)). 
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solution,  satisfy  almost  exactly  the  following  relation, 


zr9 


i9)  =  tfei9) 


sin  9 
cos  9 


(6.7) 


which  is  analogous  to  Eqn.(3.11),  as  applied  to  the  accumulated  near-tip  plastic 
strains.  Also,  as  was  observed  from  the  near-tip  strain  distribution  (Fig.  9)  of 
the  earlier  analysis,  it  is  again  found  from  the  present  computation  (Fig.  22)  that 
^g{9  =  0)  for  the  perfectly  plastic  case  is  about  3.3  times  the  corresponding  value 
given  by  the  HRR  analysis  for  large  n. 

The  near-tip  value  of  the  J  integral  was  calculated  by  substituting  epee{9)  and 
6^0 (9)  obtained  above  into  Eqn.(3.15).  The  integral  in  (3.15)  was  estimated  numer¬ 
ically,  and  it  was  found  that  Jtip  is  about  0.95  times  the  remotely  applied  J  value. 
This  is  somewhat  different  from  the  development  in  plane  strain  where  TRACEY 
(1976)  reported  Jtip  to  be  about  0.8  times  the  applied  J  value.  But  later,  SHIH 
(1981)  found  Jtip  to  be  0.96  times  the  applied  J  from  his  finite  element  calcula¬ 
tion  under  plane  strain,  small-scale  yielding  conditions  for  the  perfectly  plastic  case 
based  on  a  different  type  of  singular  element. 

If  the  near-tip  J  computed  above  from  the  present  analysis  is  used  to  normalize 
the  crack  tip  displacement  6t,  it  is  found  that  6t  =  0.88(Jtip/cr0).  Hence,  it  is 
concluded  that  <5t/(J/<70)  for  the  perfectly  plastic  case  under  plane  stress,  small- 
scale  yielding  conditions  could  vary  from  0.84-0.88. 

In  closing,  it  is  observed  that  all  the  results  given  above  by  the  present  accurate 
numerical  computation  are  in  good  agreement,  in  every  respect,  with  the  earlier 
analysis,  which  employed  a  nonfocussing  mesh  with  nonsingular  elements  near  the 
crack  tip.  The  earlier  analysis  relied  purely  on  the  fineness  of  the  mesh  and  a  large 
plastic  zone  to  the  smallest  element  size  ratio  to  provide  sufficient  resolution  near 
the  crack  tip. 


-59- 


APPENDIX 

Newton- Raphson  Method  for  Equilibrium  Iteration 

It  was  observed  in  Sec. (2)  that  an  iterative  Newton-Raphson  method  was  used 
in  the  solution  of  the  incremental  equilibrium  Equations  (2.11).  This  procedure  is 
summarized  below  for  the  kth  equilibrium  iteration  of  the  (t  4-  At)th  time  step. 

1)  The  externally  applied  load  is  increased  and  F(t  +  At)  is  calculated. 

2)  The  tangent  stiffness  matrix  K^-1(t  +  At)  and  the  vector  Pk-1(t  +  At)  = 
JrBV-1  (t  +  At)dA  are  calculated.  For  the  first  iteration  of  the  time  step 
(k=l),  the  above  vector  is  computed  from  the  converged  solution  at  the  end  of 
the  previous  time  step  as  P°(t  +  At)  =  fR  BTa(t)dA. 

3)  The  following  matrix  equation  is  solved  by  Gauss  elimination: 

K^T 1 A Uk  =  F(t  +  At)  -  Pk_1  =  ARk  . 

4)  The  nodal  displacements  and  element  strains  are  updated  as  follows, 

Uk(t  +  At)  =  Uk_1(t  +  At)  +  AUk 

ek(t  + At)  =  B  Uk(t  +  At)  . 

For  the  first  iteration  of  the  time  step  (k=l), 

U1(t  + At)  =  U(t)  +  AU1  . 

5)  In  order  to  prevent  fictitious  (numerical)  elastic  unloading  of  elements  in  some 
parts  of  the  plastic  zone  during  the  subsequent  iterations  (k  >  1)  of  the  time 
step,  a  path  independent  scheme  is  used  to  update  element  stresses.  The 
stresses  are  estimated  by  integrating  from  the  values  at  the  end  of  the  pre¬ 
vious  accepted  equilibrium  configuration  to  the  current  iteration  of  this  time 

step  by  using  the  cumulative  strains  as  follows  (BATHE  (1982)), 

rdf  (t+ At) 

ak(t  +  At)  =  £(t)  +  /  Dde 

hi  t) 
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An  explicit  method  was  employed  to  evaluate  the  integral  in  the  above  equation. 
6)  The  Euclidean  norm  of  the  out-of-balance  force  vector  ARk  (see  Step (3))  and 
the  internal  energy  increment  are  checked  for  convergence  by  comparing  with 
the  corresponding  values  at  the  start  of  the  iteration  process  as  (BATHE  and 
CIMENTO  (1980)), 

||ARk||  <  6jp||AR1|| 

AUk  .  ARk  <  ^AU1  .  AR1  , 
where  and  6e  are  small,  preset  tolerances. 

If  convergence  is  not  achieved ,  control  is  returned  to  Step(2)  to  perform  the 
next  iteration. 

If  convergence  is  achieved,  control  is  returned  to  Step(l)  to  perform  the  next 
time  step. 

Explicit  Integration  of  Incremental  Constitutive  law 

The  method  of  stress  computation  mentioned  in  Sec. (2)  is  outlined  for  an 
isotropic  hardening  solid  below. 

1)  After  solving  the  finite  element  equilibrium  equations  for  the  nodal  displace¬ 
ment  increments  AU,  the  strain  increment  Ae  is  obtained  as 

Ae  =  BAU, 

where  B  is  the  strain-displacement  matrix. 

2)  An  elastic  estimate  Acte  for  the  stress  increment  is  computed  as 

Acte  =  CAe  . 

3)  A  trial  stress  state  crE  .=  a0  +  A crE  is  calculated  from  the  stress  state  a0  at 
the  beginning  of  the  iteration.  Here  is  taken  to  be  inside  the  yield  surface 
(Fig.  23)  for  the  sake  of  definiteness. 
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Figure  23..  Stress  computation  in  the  finite  element  scheme  based  on  an  explicit 
integration  of  the  incremental  constitutive  law. 
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4)  If  F(aE)  —  (a0)2  <  0,  where  a0  is  the  value  of  a  at  the  beginning  of  the  iteration, 
then  the  elastic  behaviour  assumption  holds  and  the  remaining  steps  in  this 
method  are  omitted.  Otherwise,  the  yield  surface  has  been  crossed  during  the 
trial  stress  incrementation  (Fig.  23). 

5)  The  contact  stress  state  gc  is  obtained  as 

ctc  =  g°  +  qA<7E  , 

where  0  <  q  <  1  and  F(ar)  —  (a0)2  =  0.  This  condition  for  the  Von  Mises  yield 
function  leads  to  a  quadratic  equation  in  q.  It  should  be  observed  that  the  path 
from  a0  to  crc  constitutes  fully  elastic  material  response. 

6)  A  stress  state  aT  is  obtained  as 

2T  =  eC  +  £((l-q)A£-^F£) 

=  2C  +  ((2E-2C)-^SF£). 

In  this  equation,  F£  is  taken  as  the  normal  to  the  yield  surface  at  the  stress 
state  ac .  Also,  A  A  is  evaluated  corresponding  to  the  stress  state  gc . 

7)  The  yield  surface  is  updated  as 

oF  =  <7°  +  H(<70)  Aep  , 

where  Aep  =  |AA(7n  and  H(a°)  =  which  can  be  obtained  from  (2.4) 

for  hardening  solids  and  is  set  equal  to  zero  for  perfect  plasticity. 

8)  Due  to  the  finite  nature  of  the  time  step,  the  stress  state  aT  obtained  in  Step (6) 
will  not  (in  general)  lie  on  the  updated  yield  surface.  aT  is  then  simply  scaled 
as  follows, 


The  path  from  g}  to  <rF  constitutes  elastic-plastic  material  response. 
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In  order  to  minimise  the  error  due  to  the  use  of  finite  increments,  the  excess 
stress  £E  -  ac  is  divided  into  m  subincrements,  and  steps  (6)  to  (8)  are  carried 
out  m  times  with  the  subincrements. 
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